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Abstract. An arithmetically simple method has been developed for the evaluation of Biot- 
Savart integrals on tetrahedralized distributions of vorticity. In place of the usual approach of 
analytical formulae for the velocity induced by a linear distribution of vorticity on a tetrahedron, 
the integration is performed using Gaussian quadrature and a ray tracing technique from computer 
graphics. This eliminates completely the need for the evaluation of square roots, logarithms and 
arc tangents, and almost completely eliminates the requirement for trigonometric functions, with no 
operation more costly than a division required during the main calculation loop. An assessment of 
the algorithm's performance is presented, demonstrating its accuracy, second order convergence and 
near-linear speedup on parallel systems. 
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1. Introduction. An important part of many calculations in fluid dynamics 
and electromagnetism is the evaluation of a Biot-Savart integral, for velocity in fluid 
dynamics and for magnetic field in electromagnetism. In fluid dynamics, the source 
term is vorticity while in electromagnetism it is current. The Biot-Savart integral for 
velocity v due to a distribution of vorticity u> over a volume V is: 



where r = x — x l5 R = |r| and subscript f indicates a variable of integration. For 
electromagnetic calculations, v is the magnetic field and u> the current in the region 
V. For the purposes of this paper, it is assumed that the volume V is discretized 
into tetrahedral elements within which the source varies linearly. The question then 
is how to compute the resulting field v. This paper is motivated by the requirement 
to evaluate velocities in Lagrangian vortex methods, where a moving distribution of 
control points is tetrahedralized at each time step as an aid to velocity calculation. 
Such a method has been implemented by Marshall et al. [4], using a mixture of ana- 
lytical formulae and Gaussian quadrature to carry out the integration of equation 1.1. 
The aim of this paper is to develop a method which is arithmetically simpler than 
that used hereto. 

A number of analytical formulae and numerical procedures have been developed 
for the evaluation of equation 1.1. A recent general paper is that of Suh [11] who 
applies Stokes' theorem to reduce the volume integral over an element to a number 
of surface integrals which are in turn reduced to line integrals which can be evaluated 
analytically. This work is similar to that of Newman [6] who derived equivalent 
formulae, focussing on applications in fluid dynamics. In electromagnetism, there is 
an extensive literature on the evaluation of equation 1.1, including analytical [2, 7, 12- 
14] and numerical [3] approaches. As this paper is motivated by the fluid dynamical 

"This work was supported by the European Community-funded project HPC-Europa, contract 
number 506079, and was carried out at the Institut fur Stromungsmechanik und Hydraulische 
Stromungsmaschine and HLRS, both at the University of Stuttgart, Germany 

t Department of Mechanical Engineering, University of Bath, Bath BA2 7AY, United Kingdom 
(m . j . carleyObath . ac . uk) 




(1.1) 



2 



M. CARLEY 



problem, it will use the terms 'velocity' and 'vorticity' but it should be noted that 
there are also many useful related results in the electromagnctism literature. 

In evaluating the velocity at the nodes of a vorticity distribution which is dis- 
cretized into tctrahcdra, a Gaussian quadrature can be used for tetrahedra which 
are far from the evaluation point x, but (integrable) singularities in the integrand 
make this awkward for nearby elements. In this case, the standard approach is to 
use analytical formulae which have the benefit of being exact and non-singular. The 
disadvantage, as noted by Marshall et al. [4], is that such formulae [6] are compu- 
tationally expensive, requiring in this case 12 logarithms and 24 arc tangents per 
tetrahedron. Even if the analytical formulae are only used for 'nearby' elements, they 
still represent a large part of the velocity computation. The aim of this paper is to 
present a velocity computation method which is arithmetically simple, requiring no 
operation more onerous than a division during the main calculation, by using a ray- 
tracing method borrowed from computer graphics. The resulting method can then be 
used in codes as a 'plug-in' replacement for previous techniques. 

2. Velocity evaluation. The integral of equation 1.1 is to be evaluated at a 
number of points x, which also form the nodes of a distribution of vorticity u>. In this 
case, as in the work of Marshall ct al. [4], the nodes are tetrahedralized so that they 
form the vertices of a collection of tetrahedra. It is further assumed that vorticity 
varies linearly over the elements. The aim is now to avoid the computational effort 
involved in the standard analytical approach to quadrature and develop a method 
which is arithmetically as simple as possible. In the words of Richardson, this is 
a problem where it may be quicker to arrive at a destination at the "footpace of 
arithmetic" , rather than "on the swift steed of symbolic analysis" [8] . 

First, equation 1.1 is rewritten in spherical polar coordinates centred on the eval- 
uation point x: 

v=— / / / s x wdRsin^d^dfl, (2.1) 
4tt J J J_ 00 

where r = x — xi = — Rs with the unit vector s = (sin ^cos 9, sin 0sin#, cos <fi). This 
transformation also eliminates the integrable singularity in the integrand, making the 
integration rather easier from a numerical point of view. The azimuthal and polar 
integrals are evaluated using Gaussian quadratures so that: 

N M w^wHt f°° 

v » ^2 X! sin ^" "i a T \ SnmXudR, (2.2) 

n=lm=l J -°° 

where 

&, = (l+tW)7T/2, 
ro = (1 + 4^/2, 

§„m = (sin <f>„ cos 9 m , sin <j> n sin 6 m , cos <j> n ) , 



and {t\" , ,w ( i > ), i = l,...,N, are the nodes and weights of an N point Gaussian 
quadrature with —1 < t < 1. The main computation is now reduced to integrals over 
R along rays in the direction s. 
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2.1. Integration over tetrahedra. The radial integral of equation 2.2 must 
now be evaluated, to include the contribution of the mesh tetrahedra. The approach is 
summarized in figure 2.1: if the ray from x in the direction s intersects the tetrahedron, 
i.e. passes through two of its faces, the tetrahedron makes a finite contribution to the 
integral. The overall integral is: 

/oo 
sxwdi?. (2.3) 
-oo 

If the ray enters a tetrahedron at R = R n and exits at R = R\, 

I{x,6,4>) = Y J [ Rl s*"dR =i^( J R 1 -i? )§x (« +wi), (2.4) 

J Rq 

summing over tetrahedra which are intersected by the ray. The vorticities u> and u>i 
are those at the entry and exit points Rq and R\ and, as noted previously, oj varies 
linearly between them, so that equation 2.4 is exact. If the ray enters or exits a face 
of the tetrahedron at area coordinates (u, v), the vorticity u) is: 

UJ = (1 — U — 1))Wq + MWl + VLJ 2 , 

where u>i 1 i = 0,1,2 are the vorticities at the triangle vertices. This integral is 
arithmetically very simple, with no operation more complex than a division, not even 
requiring a square root for the calculation of a distance. 

The remaining part of the algorithm is how the intersection, if any, of a ray 
with a tetrahedron can be determined. This is found using a method from computer 
graphics. Moller and Trumbore [5] give a method and C source code for determining 
the area coordinates (u, v) and directed distance R of the intersection point of a ray 
with a triangle. Their method is very efficient requiring at most one division and no 
operation more complex than this. It is used in the algorithm of this paper to find 
the intersections of a ray with a tetrahedron, by checking the faces in turn. In order 
to accelerate the procedure, an initial check is carried out on the tetrahedron as a 
whole. 
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Fig. 2.2. Pre-check to decide if tetrahedron is likely to be intersected by a ray 

Figure 2.2 shows the check for the case of a triangle. The vector r is the vector 
from x to an arbitrary vertex of the tetrahedron. The scalar product of s and r is: 

r.s = |r| cosd, 

where 8 is the angle between s and r. If the tetrahedron subtends less than this angle, 
none of its faces are intersected and no further check is needed. If it has a maximum 
edge length h, the maximum angle it can subtend is ip w h/\r\ and cos 6* < cosi/> is 
a necessary condition for the tetrahedron to be intersected. Using the small angle 
approximation of cos?/; and neglecting fourth order terms: 

cos 8 < 1 - h 2 /\r\ 2 , 
cos 2 6< 1 -2h 4 /\r\\ 

so that the condition can be written: 

(r.s) 2 < \r\ 2 -2h 2 . 

The test is implemented by setting h 2 to the square of the longest edge length on the 
element and by applying it only for |r| greater than some minimum value. 

2.2. Summary of algorithm. In summary, given a set of points Xj each with 
an associated vorticity u>j, the Biot-Savart integral at each point can be evaluated as 
follows: 

1. Generate a tetrahedralization of the points Xj using, for example, the TET- 
GEN code of Si [10]. 

2. Select Gaussian quadrature rules of order TV and M for integration in <f> and 
8 respectively. 

3. Precompute the factors k nm — sin (f> n w n N ' 'wm^ /I6ir 2 and the ray directions 
s nm = (sin 4> n cos8 m , sin <f> n sin# m , cos <p n ), equation 2.2, for n = 1, . . . , N and 
m=l,...,M. 

4. For each tetrahedron, loop over the nodes Xj and check for intersection of 
each ray s in turn. If the ray intersects the tetrahedron, add the contribution 
of equation 2.4 to the velocity Vj. 
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Table 3.1 

Computation time and r.m.s. error in velocity calculations for Hill's spherical vortex 



N 


T/s 


e 


T/s 


e 


T/s 
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1000 


2 


0.200 


28 


0.036 


13 


0.036 


2000 


8 


0.100 


100 


0.028 


44 


0.026 


4000 


30 


0.020 


369 


0.021 


155 


0.017 


8000 


100 


0.006 


1385 


0.019 


562 


0.013 


16000 


400 


0.004 






2168 


0.009 




Cray C90 


Laptop 


8-nodc cluster 



If the calculation is being carried out on a parallel system, the tetrahedra can be 
divided amongst the processors for the velocity calculation with the total velocity 
being found by a summation at the end of the computation. This only requires one 
communication per velocity calculation. 

3. Numerical tests. A number of tests have been carried out to evaluate the 
performance of the method of §2.2 with respect to speed, accuracy and paralleliza- 
tion. The method has been implemented in a Lagrangian vortex code currently under 
development, similar to that of Marshall et al. [4]. The control points are tetrahedral- 
ized using the TETGEN code of Si [10] and the program runs on serial and parallel 
systems using the MPI message passing interface. 

3.1. Speed and accuracy. The first test case considered is Hill's spherical 
vortex [9, pp23-25]. This is a distribution of vorticity with a linear variation of 
azimuthal vorticity ujg = Ar inside a sphere of radius a. Within the sphere, radial 
and axial velocities are: 



This is an especially useful test case because, with the exception of errors in dis- 
cretizing the boundary of the sphere, the vorticity distribution is exactly duplicated 
by linear interpolation over elements. The test was carried out by evaluating the 
velocity at N randomly placed points with a sphere of radius a = 1 with N varying 
from 1000 to 16000. Calculations were carried out a 500MHz Intel Pentium III laptop 
with 4 point Gaussian quadrature in <f> and 6 and with 16 point quadratures on an 8 
node 3.2GHz Intel Xeon cluster. Table 3.1 shows the r.m.s error in velocity and the 
computation time for these two cases and, for comparison, the corresponding data for 
Marshall et al.'s calculations on a Cray C-90 [4]. From the data presented, it appears 
that for small numbers of nodes, the method is more accurate than the use of analyti- 
cal formulae and that it remains better up to distributions, in this case, of about 4000 
points. The second point to note is that the accuracy appears to be controlled, up to 
the limit of the resolution of the vorticity, by the number of quadrature points. The 
4x4 quadrature error stops reducing at about 4000 nodes while the 16 x 16 is still 
improving at 16000 nodes. 

3.2. Convergence. As a second test on the accuracy of the calculation method, 
and to assess the convergence properties, the velocity due to a Gaussian core ring was 
computed. This has exponentially small vorticity on the boundary of the vorticity 
distribution so that errors due to discretization of the boundary should be reduced, 
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Fig. 3.1. Convergence of velocity for a Gaussian vortex ring, r.m.s. error against number of 
quadrature points: solid line 16 X 17 points; dashed line 64 X 65 points; dotted line 128 X 257 points. 




0.025 



Fig. 3.2. Error in velocity against discretization length h, for 64 X 64 point quadrature: circles: 
errors; solid line: linear fit e = 8.03/i 2 15 



compared to the case of Hill's spherical vortex. For comparison, the velocity was 
computed using the stream function for an axisymmetric vortex ring [9, pages 192- 
194], differentiating it numerically on a dense grid to give the velocity. This velocity 
was interpolated to find the velocity at the control points of the test distribution. 
Three different distributions of points were considered, each made up of a number of 
stations equally spaced in azimuth, with a square grid of points at each station. The 
low resolution distribution had 16 x 17 points, the intermediate 64 x 65 and the high 
resolution 128 x 257. The velocity was evaluated using equal numbers of quadrature 
points in cj> and 9. Figure 3.1 shows the r.m.s. error in velocity as a function of the 
number of quadrature points for all three test cases. 

It is clear that, as proposed in §3.1, the accuracy of the integration is controlled 
by the order of Gaussian quadrature, up to a limit fixed by the resolution of the 
vorticity distribution. This is investigated in figure 3.2 which shows the r.m.s error 
e of figure 3.1 plotted against a typical discretization length scale h. The velocities 
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N 

Fig. 3.3. Computation time T for 16000 point Hill's spherical vortex as a function of number 
of processors N: circles computation time, solid line linear fit. 

were computed using 64 quadrature points in <fi and 9. At this stage the two lower 
resolution test cases have reached their minimum error. The error scales as h 2A5 
showing that the method is approximately second order. 

3.3. Parallelization. The final requirement of a velocity calculation algorithm 
is that it be possible to implement it efficiently on a parallel system. The method of §2 
has been coded for an MPI system by sharing the tetrahedra out amongst the proces- 
sors and combining the velocity contributions at the end of the calculation. Only one 
communication is needed, to combine the velocities computed on each processor so 
that the network overhead is low. Figure 3.3 shows the computation time for velocity 
on Hill's spherical vortex discretized with 16000 points, with 4-64 3.2GHz processors. 
The fitted line shows the calculation time scaling as jV -0 - 98 , a near linear speedup. 

4. Conclusions. A method for the evaluation of Biot-Savart integrals has been 
presented which eliminates the need for complicated mathematical operations over 
most of the calculation. With the exception of a negligible number of trigonometric 
factors which must be pre-computed, the method requires no operation more costly 
than a division. Calculations performed using the technique demonstrate that it is 
approximately second order accurate, that it converges correctly and that it achieves 
near-linear speedup on a parallel system. Future work will consider how to use the 
approach in methods similar to fast multipole and how to improve the ray trac- 
ing algorithm used, perhaps by taking advantage of adjacency relationships between 
tetrahedra. 

In conclusion, we note that the method is based on a ray-tracing technique from 
computer graphics and that such methods are now being implemented at hardware 
level on mass market systems. This opens the possibility of using the calculation 
method on relatively cheap hardware, giving further acceleration at little additional 
cost. 
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